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Abstract 

A phenomenological model based on the three-dimensional theory of 
nonlinear elasticity is developed to describe the phenomenon of over- 
stretching in the force-extension curve for dsDNA. By using the concept 
of a material with multiple reference configurations a single formula is 
obtained to fit the force-extension curve. 

Keywords: force-extension curve, overstretching, nonlinear elasticityjimiting chain 
extensibility 

Abbreviations: dsDNA, double stranded DNA; ssDNA, single stranded DNA; 
WLC, Worm like chain. 

1 Introduction 

A typical force-extension curve for dsDNA exhibits three portions [H [H [3] . Dur- 
ing the first portion there is an entropic stretching regime (usually modeled by 
the worm-like chain), followed by a force plateau in the region of 65 picoNewtons, 
while in the last portion there is a sharp transition from the usual B-form to 
a new overstretched form, usually designated S-DNA. The structure of S-DNA 
remains the subject of debate but it should not be confused with ssDNA [H[S]. 



The biological function of the overstretching DNA transition is complementary 
to thermal or pH induced denaturation and for this reason there is considerable 
attention focused on modeling this phenomenon [5] , but there appears to be no 
general agreement about the models that have been proposed in the literature. 

In [7] a model based on a force-induced melting of the DNA double he- 
lix was proposed. This implies that S-DNA is made up of a mixture of large 
islands of separated ssDNA and remnant base-paired B-DNA, and molecular 
extension is a weighted average of its extension in the two possible states. A 
two-state worm-like chain has been also proposed in [5]. In several papers the 
idea that the overstretching behavior of DNA may be modeled by a sort of 
mixture theory has been applied in the study of B-DNA to S-DNA transition 
as a function of solution conditions, including variations in temperature, pH 
and ionic strength (see, for example, [9]). In [10] . by using a thermodynamical 
model for tension-melted dsDNA it is argued that the overstretching transition 
cannot be explained in terms of conversion of double helix to nonintcracting 
polynucleotide strands. This is because two parallel noninteracting ssDNAs 
cannot explain quantitatively the mechanical properties of S-DNA. This is ar- 
gued directly from an examination of the experimental data by the authors [10] . 
The Rouzina and Bloomfield model [7J [5] is therefore criticized because in the 
B-ss scenario the overstretched state should be associated with a constant force 
between the B-DNA and ssDNA lower than that observed. 

The aim of the present note is to present a new framework for describing 
the overstretching phenomena. Our model is developed using a non-standard 
version of the phenomenological theory of nonlinear elasticity where the stress 
is determined as a function of the deformation gradient calculated with respect 
to a varying reference configuration in a such a way that it is possible to in- 
troduce micro-mechanical considerations. This idea was introduced originally 
by Eckart and then developed more recently by Rajagopal and Wineman 
[12l [13] in order to formulate constitutive equations for materials that undergo 
deformations induced by microstructural changes. Recently De Tommasi et al. 
[14] have proposed a micro-mechanical interpretation of this theory that may 
be quite useful in the study of the overstretching phenomenon. 

The general form of the constitutive equation for nonlinear elasticity is ex- 
pressed in terms of a strain-energy function. In the standard theory it is as- 
sumed implicitly that the material response is due to a molecular mechanism 
that does not change during the deformation process under consideration. In 
single-molecule experiments on DNA, this assumption may be considered valid 
only on the first portion of the force-extension curve. At a certain moment 
the hydrogen bonds between strands start to break and there is a fundamental 
change in the molecular mechanism responsible of the overall material response. 
In the case of DNA these microstructural changes are driven by several factors: 
stretching, salinity, temperature, etc. From a micro-mechanical point of view 
it is possible to look at the nucleotides as particles that are connected by two 
different types of chains. A fraction of the chains is elastic and endows the 
DNA molecule with nonzero stiffness. The complementary chains are break- 
able and are responsible for the alteration of the molecule. The stress in each 
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breakable chain is zero until a certain activation threshold is reached and af- 
ter a limiting value of the strain is overcome. We assume that a continuous 
process of microstructural conversion occurs after the deformation increases be- 
yond a threshold value. In this initial model we shall neglect any factors such 
as salinity, temperature or ionic strength that are not strictly mechanical. We 
emphasize that instead of considering the conversion of the double helix to 
two noninteracting polynucleotide strands, we are considering here a process of 
conversion related to the rupture of stress-bearing bonds. Upon rupture of the 
bonds a new microstructural arrangement forms with a new unstressed reference 
configuration. More details of this constitutive model may be found in [13] . 

Let us consider a deformation x = x(X,i), where x is the current position 
of a particle located at X in the undeformed configuration at time t = 0. The 
deformation gradient is given by F(X, t) = <9x/9X and the left Cauchy-Green 
tensor by B = FF T . We assume that there is a range of deformation for which 
the material behaves like an incompressible, isotropic elastic material, i.e. the 
Cauchy stress T = —pi + T*^' 1 ), where —pi is the indeterminate part of the 
stress due to the constraint of incompressibility (detF = 1) and the extra stress 
takes the form 

T^' 1 ) = 2IU 1 (1) B - 2T^ 2 (1) B- 1 . (1) 

The strain-energy function = W^(Ix,I 2 ) is a function of the principal 

invariants I x = tr(B) and I 2 = tr(B- 1 ), where wf 3 = dW^/dI i: i = 1,2. 
An activation criterion is needed to determine when the microstructural change 
begins. This is provided by introducing a scalar deformation state parameter s. 
Here, we suppose that s = s(I±, I2) depends on the deformation through I\ and 
I2 for consistency with the requirement of isotropy, although more general forms 
of s may easily be adopted. For s < s a , the threshold value of s, no conversion 
has yet occurred, i.e. all the material is in its original form and the stress is 
given by {1}. On the other hand, for a value of the state parameter s~ beyond 
s a microstructural changes have occurred and the reference configuration has 
changed. This implies that the stress is now a function of the relative defor- 
mation gradient for the material formed at state s given by F = <9x/<9x, where 
x is the position of the particle in the configuration corresponding to deforma- 
tion state s~. In Figure 1 the original reference configuration, the configuration 
at J and the current configuration are depicted. The associated Cauchy-Green 
tensor is given by B = FF T . 

We shall assume that the new material formed at the state s is still elastic, 
isotropic and incompressible such that the extra Cauchy stress at state s in this 
new configuration formed at the deformation state s is given by 

T (E,2) _ 2W[ 2) % - 2W r 2 (2) ET 1 . (2) 

Here = W^\h,I 2 ) is the strain -energy function of the newly formed 

material, relative to the reference configuration at 5T. Another important sim- 
plifying assumption is that a single function governs the strain energy 
during the continuous microstructural change. The total current stress is taken 
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Figure 1: Schematic of a material with an evolving reference configuration 

as the superposition of the contributions from the material remaining in its orig- 
inal configuration and from all the new material formed at deformation states 
s e [s a ,s], i.e. 

T=-pI + b(s)T^ + J a(s)T {E ' 2) ds. (3) 

In ([3]) the function a(s) is a conversion rate satisfying a(s) = when s < s a 
and a(s) > for s > s a , while b(s) is the volume fraction of the material in 
the original configuration remaining at state s, with 6(s) = 1 when s < s a and 
< b(s) < 1 for s > s a . Thus, to complete the model, constitutive equations 
for and W^ 2 \ the activation criterion and the conversion rate have to be 

prescribed. Our model is three-dimensional and fully consistent with the theory 
of continuum mechanics. To illustrate the ideas quantitatively we begin with 
a prototype that is empirical and one-dimensional, and we then show how to 
recast the theory in three-dimensional form. 

2 The constitutive models 

2.1 Data sources 

We consider the sets of experimental data in [H [5] , which correspond to dif- 
ferent salt concentrations. We use measured data (xi,fi), i — l,...,m, here 
corresponding to a force- (/ in picoNewtons)-extension (x in microns fx) exper- 
iment on a single dsDNA molecule in 250 mM [Na + ] buffer solution at 7.5 pH 
(these data are reported in figure 3 of [4]). 
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2.2 An empirical one-dimensional model 

Let x denote the one-dimensional extension. On the same basis as illustrated 
above, for the one-dimensional force / we have 

f{x)=b(x)f<»+ I* a(x)fW(x)dx. (4) 



In the various quantities a(x),b(x),x a have the same meaning as before, 
with x replacing s. Since the process of conversion is continuous we have 

b(x) = 1—1 a(x)dx, x > x a . (5) 

The constitutive assumptions we introduce are: for a logistic modifi- 
cation of the original one-dimensional Fung model widely used in biomechanics 

ma, i.e. 

f (D M _Mi exp[/3(a - sp)] 

1 [X) 2 exp[/?(x-*o)]+ 7 ' [ > 

where the material constants Hi, Xg and [3 have dimensions of force, length and 
1/length, respectively, and 7 > is a dimensionless constant; and, for f^ 2 \ the 
WLC interpolation formula 



f {2) (x)=V2 



1 / x-2 1 



(7) 



Here ^2 = kT/l p , where l p is the persistence length, k is Boltzmann's constant 
and T the temperature (degrees Kelvin), z — x/l c and l c is the contour length of 
the molecule. We have chosen a logistic Fung model to capture the first portion 
of the force-extension curve (i.e. to capture the strain-hardening phenomenon) 
but without introducing a singularity such as that in the WLC formula. The 
WLC is used to model the sharp increase in force at the end of the curve just 
after the plateau. It is clear that the modeling of the plateau zone depends 
on how the reference configuration evolves, and this may be controlled by the 
choice of the conversion rate. Usually, in the context of rubber mechanics, very 
simple models for the conversion function are adopted (for example, quadratic 
or piecewise linear functional forms). Here we use a functional form suggested 
by statistical mechanics, namely a probability distribution function computed 
by considering two possible states for a chain composed of a fixed number of 
base pairs with a given fixed difference in the energy between the two states. 
For this purpose let 

9 ^ X ' ~ [l+e-ci(z-c 2 )]2' 

where ci, C2 and 5 are constants, and define 

a(x) =g(x) -g(x a ) xe[x a ,x c ], (8) 
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with a(x) — otherwise. Here, we are assuming that the conversion has been 
completed when x reaches the value x c , and this imposes the continuity require- 
ment a(x c ) = 0, which leads to C2 = [x a + x c )/2. A plot of the function a(x) is 
shown in Figure 2. 

We denote by C the total fraction of the material that can undergo conver- 
sion. Then, 



C= a(x)dx. (9) 
From the definition of C in ©, we calculate 

[ 1 + e -c 1 (x„- Ca )][ 1 + e -c 1 (x.-c a )] 

O/C = 7 ; -, ; q(x a )(X c — X a ). 

1 g-ci (x a ~c 2 ) _ g-C! (x c -c 2 ) yv a;v c ay 

The constitutive parameters to be found in this empirical model are /ii, /?, 7, /i2 
and Z c . Moreover, we have to fix the activation criterion and therefore we also 
need values for x a , x c , c\ and C. At this stage the only a priori information 
about these parameters is that C £ [0,1]. The strategy for fitting that we use 
to deal with the original force-extension data (X4, fi),i = 1, . . . , m, is explained 
in the Appendix. Since, in principle, several parameters have to be identified 
in our model, their numerical approximation could pose severe problems (see, 
e.g., |16|). For this reason we devise a strategy that accounts for the physical 
interpretation of some of these parameters. A set of parameters identified by 
the fitting results is given by 

x a = 20, x c = 28.754, C = 0.50367, 

m= 64.977, /3 = 2.7537, 7 = 0.019288, } CIO) 
ci = 0.059145, /i 2 = 25.87, l c = 32.022, 

with residual res 2 = 21.93. 

In Figure 3 the prediction of the model obtained by using these parameters 
is shown. The results are quite good, but we believe that better insight might 
be gained from the three-dimensional model. 

2.3 Three-dimensional models 

In three dimensions the single molecule force-extension experiment is idealized 
as a simple tension test, for which the deformation is given by 

where A is the stretch in the axial (i.e. z) direction. The current deforma- 
tion gradient is given in matrix form by F(A) = diag(l/\/A, 1/V% A), and the 
corresponding Cauchy-Green deformation matrix is B(A) = diag(l/A, 1/A, A 2 ). 
Hence, 

/i=A 2 + 2A-\ J 2 =A" 2 + 2A. (12) 
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Figure 3: Plot of f(x) vs. x. Circles are from the data in Wenner et al. |4J 
for 250 mM. The curve is obtained from the fitting procedure, which yields the 
parameter values given in (|10|) with squared residual res 2 — 21.93. The red stars 
indicate the range of values for which the conversion is active 
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Since this is a one-parameter deformation, it is possible to establish that 
there is a one-to-one correspondence between the activation parameter s and the 
stretch, and we write s = s(X). For this reason we use the terminology activation 
stretch, which we denote by A a , instead of a generic activation parameter s a 
(then, s a = s(X a ) — X a ). The deformation gradient at state A is therefore 
denoted by F(A) = F(A)F _1 (A) and we therefore compute 



F(A) = diag I y A/A, y A/A, A/A 

B(A) = diag (A/A, A/A, A 2 /A 2 ) . (13) 



It follows that, for example, I\ = A 2 /A 2 + 2A/A. If we consider the class of elastic 
materials referred to as generalized neo-Hookean materials, with W = W(Ii), 
then from (fTJ) we obtain the principal components of the Cauchy stress tensor 
in the form 

U = 2X^W 1 -p, i = 1,2,3. (14) 

The requirement that the lateral surfaces of the specimen undergoing the simple 
extension are traction free, t± = ti = 0, yields 

V = 2X- 1 W 1 . (15) 

Generalizing these results to the case ^ the tensile force per unit deformed 
cross-sectional area necessary to achieve the stretch is given by the Cauchy 
stress component 



t 3 (A) = 26(A) (A 2 — A -1 ) W{ 



(i) 




w[ 2) {X)dX, (16) 



where 



A 



6(A) = 1 - / a(X)dX. (17) 



The corresponding force per unit undeformed area of cross-section is A -1 £3 (A). 

At this point it is necessary to complement (|16|) with the constitutive equa- 
tions. We need a constitutive equation for the strain-energy function of the 
material before the conversion starts, i.e. W^ 1 **, to model the first portion of the 
force-extension curve. Then, we also need a constitutive equation for the func- 
tion that governs the mechanical behavior of the newly formed material. 
This choice is important for modeling the "last" portion of the force-extension 
curve. The overstretching plateau, as already pointed out, is modeled by the 
choice of the conversion function a(s). For the strain energy in the first 
regime we consider a modification of the strain-energy function, here denoted 
W F , proposed by Fung for modeling biological tissues. This is given by 

= |«pt9(J 1 -3)]. (18) 
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As for the ID case, we need to modify this relationship because a saturation 
phenomenon has to be taken into account. The mechanical behavior charac- 
terizing the strain stiffening of the DNA molecule in the first portion of the 
force-extension curve cannot influence what happens in the plateau zone. For 
this reason we consider a logistic modification of the (three-dimensional) Fung 
model (fTgj) analogous to that used for ID. This is given by 

w {i) _ Mi expWi ~ 3 )] / 1q x 
2 exp[/3(/ 1 -3)]+ 7 1 

so that 

^ (1) = ^ln(exp(/3(/ 1 -3))+ 7 ), (20) 

which reduces to the neo-Hookean material W — ~ 3)/2 when 7 = 0. For 

the strain-energy function in the second portion of the deformation range 
we consider the phenomenological model first proposed by Gent [17] and given 

by 

W (2) (/i) = -^J m ln^l-^_y h<J m + 3, (21) 

where \i% is the shear modulus for infinitesimal deformations and J m (> 0) is 
the limiting value of I\ — 3 associated with limiting chain extensibility. In the 
limit as the chain extensibility parameter tends to infinity ( J m — > oo), (I21|) also 
reduces to the classical neo-Hookean model. The model (|2H has been discussed 
in detail by Horgan & Saccomandi [TB] and it can be connected with the so- 
called Freely Jointed Chain (FJC) model. In this case the response function is 
given by 

(2) = /£ Jm (22) 

2 j ro _ (7l _ 3) ' 

so that the stress has a singularity as I\ — > J m + 3. 

The model we have proposed contains several constitutive parameters that 
have to be found by using a fitting procedure. The parameters needed to fix 
the strain-energy functions are 7 and /i2, J m - Moreover, we use the same 

activation criterion as was used for the ID model in (JSj> . Hence, to fix the 
activation criterion we need to identify the interval with ends x a — > A and x c — > 
A c and the parameter c\ . Note that this criterion may easily be reformulated in 
a way compatible with 3D elasticity in terms of the invariant 1\ . Equation (j 16[) 
provides a formula for the Cauchy stress, but it is the nominal stress \~ 1 t 3 (\) 
(force per unit reference cross-sectional area) that is needed for the data fitting. 
We therefore transform the data set (xj,/,) into the data set (Aj,/j), where 
Ai = 1 + Xi/l c , with the contour length l c identified in the ID case. To match 
the dimensions of the force / in the data the stress A -1 £3 (A) has to be multiplied 
by the reference cross-sectional area, which is unknown. However, this is just a 
multiplicative factor that is accounted for by incorporating it into the constants 
fii and fi2i which then have dimensions of force as in the ID situation. 
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The parameters obtained by the strategy explained in the Appendix are 

A a = 1.3833, A c = 1.895, C = 0.71658, ~\ 

Mi =60, /3 = 43.772, 7 = 1.193e5, I (23) 

ci = 9.6293, ii 2 = 5, J m = 1.8819, J 

with squared residual res 2 = 57.68. The result of this fitting is shown in Figure 
4. 
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Figure 4: Plot of A _1 ts(A) vs. A. Circles are from the data in Wenner et al. 
[3] for 250 mM. The curve is obtained from the fitting procedure, which yields 
the parameter values given in (f2"3"| with squared residual res 2 = 57.68. The red 
stars indicate the range of values for which the conversion is active. Note that 
the dimensions of the stress have been converted into 'force' by multiplication 
by the unknown reference cross-sectional area by incorporating this into the 
parameters \x\ and n 2 

The model proposed herein gives good results in fitting the data, and be- 
cause it has been formulated within a very general framework it may easily be 
extended to take into account several variables of biological interest. The model 
is interesting not only because it is comprises a single formula describing the 
complete force-extension curve, but also from a conceptual point of view. In- 
deed, as has been argued by Cocco et al. |10) . S-DNA cannot be described as 



10 



a simple sort of mixture between the dsDNA and ssDNA. The relationship is 
more complex and it is clarified by the existence of multiple references config- 
urations. We point out that because DNA overwinds when it is stretched, we 
need a three dimensional model to obtain a complete and realistic picture of the 
single molecule experiments and our model is just a rigorous version of the toy 
model proposed by Gore et al. in [TO] . 

3 Appendix 

The strategy for fitting the theoretical model to the experimental data is based 
on a nonlinear least squares (LS) approximation as follows. As a first step, we 
fix a priori some parameters from simple biological considerations and we solve 
the optimization problem for the remaining parameters in order to identify a 
first optimal subset, p* say. In the successive steps, the strategy consists of 
implementing the LS algorithm by starting from this solution and then moving 
in a descent direction by including each time a new free parameter from amongst 
those that were fixed. The solution found at each step is then used as an initial 
guess for solving the next LS problem in which a further parameter has to be 
identified. Only x c (for the ID model) and A c (for the 3D model) are always 
fixed. This assumption implies that in the final part of the experimental curve 
for x > x c (or A > A c ) the material is all converted to its new form. 

All the computations are performed in Matlab with the Isqcurvefit routine 
(see [20]) for solving nonlinear least squares problems. We allow the algorithm 
to perform a maximum of 3000 iterations and stop with stringent tolerances on 
the errors (tol= le— 12). 

For the ID model, in the first step we use the optimization procedure to 
identify the parameters p = [fix, (3, "f) T , while the others are fixed by considering 
the following physical features: 

- the total contour length l c is chosen to be slightly larger than the last datum 
value for the extension since its value locates the asymptote of the WLC; 

- xq = x a since the value of xq in the logistic function ((6]) identifies the point 
where the largest growth occurs, and this corresponds to the meaning of x a in 
the activation criterion; 

- since fix corresponds to the horizontal asymptote of the logistic function, its 
starting guess is set to almost the force value of the plateau in the data; thus, 
we set fii = 69; 

- we set x a = 20 so that up to the beginning of the plateau the material is all 
in its original form. Moreover, we set C — 0.5, requiring by this assumption a 
conversion of 50%. 

The optimal parameter set identified in this first step, is p* = \p*, /?*, 7*] T = 
[69.0883, 2.907, 0.0153] T . Hence, by using this first approximation, we define a 
new sequence of optimization problems, where the fixed parameters are consid- 
ered free in the (arbitrary) sequence [C, X a , l c , ci, /J®]. The optimal final result 
is reported in the text and in Figure 3. 

The same fitting strategy used for fitting the data with the ID empirical 
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model (U]) is used for the 3D model (fll)]) . For the activation criterion we fix 
A c = 1.895. Moreover, we set C = 0.8. For the Gent material we fix \ii 
equal to almost the force corresponding to the plateau, and the parameter J m , 
accounting for the asymptote location, such that J m is almost the last numerical 
value available for the stretch data. At the first stage of the fitting the free 
parameters are again those of the Fung model and if p* is the set identified in 
this step, the (arbitrary) sequence in which the other parameters are considered 
as free is [C, J TO , A a , ci, /X2]. In this way a better (lower residual) optimal solution 
is found and the result is reported in the text and in Figure 4. 
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